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Driven optical lattices permit the engineering of effective dynamics with well-controllable tun¬ 
neling properties. We describe the realization of a tunable a Chern insulator by driving particles 
on a shaken hexagonal lattice with optimally designed polychromatic driving forces. Its implemen¬ 
tation does not reqnire shallow lattices, which favors the study of strongly-correlated phases with 
non-trivial topology. 


The manipulation of quantum systems has reached lev¬ 
els of accuracy that allow controlled variation of their 
properties. This permits extremely precise investigations 
of physical processes and opens up the opportunity to de¬ 
sign materials for scientihc and economic applications. In 
particular, the last years have witnessed an increasing in¬ 
terest in topological phases of matter largely motivated 
by promising applications such as topological quantum 
computing [1] or physical phenomena such as the quan¬ 
tum spin Hall effect [2-4]. 

Systems with topological properties can be classified 
with topological invariants whose non-zero values indi¬ 
cate non-trivial phases [5] . Pioneer work in this field was 
done in the context of the integer quantum Hall effect, 
when it was shown that the quantization of the Hall con¬ 
ductance is directly related to the first Chern invariant 
[6]. Thereafter, Haldane demonstrated with a model of 
particles tunneling on an hexagonal lattice [7] that a net 
magnetic field is not necessary for the quantization of 
the Hall conductance, owing to the topological nature of 
the model. Remarkably, the Haldane model has been re¬ 
cently experimentally implemented with shaken optical 
lattices [8]. Nonetheless, the exploration of its topologi¬ 
cal diagram crucially requires shallow lattices since pre¬ 
existing next-nearest-neighbor tunneling in the undriven 
system is necessary. 

Periodically driven lattices offer an extraordinary plat¬ 
form to engineer controlled dynamics. A number of the¬ 
oretical [9-15] and experimental [8, 16-19] works have 
demonstrated the possibility to modify the system dy¬ 
namics in a controlled fashion and to generate diverse ef¬ 
fects that include coherent destruction of tunneling [16] 
and the creation of synthetic magnetic fields [17-19] or 
topological properties [8]. 

Even though the effective dynamics that a driven sys¬ 
tem undergoes crucially depend on the specific time- 
dependent driving force, rather simple driving forces are 
usually employed. Nevertheless, as demonstrated in var¬ 
ious fields, including chemistry [20-22], nuclear magnetic 
resonance [23, 24], quantum information [25, 26], and 
many-body systems [27], essentially any desired dynam¬ 
ics can be induced with the appropriate choice of poly¬ 
chromatic driving at desired instances of time or during 
an extended time-window [28, 29]. 

In this Letter we show an optimal implementation of a 
Chern insulator by driving non-interacting particles on 


an hexagonal lattice with polychromatic driving. We 
demonstrate how the entire topological diagram can be 
accessed with deep optical lattices and suitably chosen 
driving forces. 

Driven systems can be used as quantum simulators due 
to the possibility to approximate their dynamics in terms 
of a time-independent Hamiltonian. According to Flo- 
quet theorem [30], the time-evolution operator of a peri¬ 
odic Hamiltonian H{t) = H{t + T) can be written as 

U{t) = 4(t)e-*^=«*17F(0), (1) 

where Upit) is a T-periodic unitary and iJeff defines a 
time-independent effective Hamiltonian. The distance 
between the exact dynamics induced by H(t) and the 
approximate dynamics C4ff(t) = induced by iLeff 

is bounded llD(t)-C/eff(t) 11 < lll-t/;(f)ll + lll-C/F(0)ll. 
Consequently, if the unitary Upit) is sufficiently close to 
the identity during an entire period, the dynamics of the 
system can be very well captured by UeS for all times. 
This condition is typically satisfied in a suitable fast¬ 
driving regime, where the driving frequency lo = 27r/T 
is the largest energy scale of the system. The effective 
Hamiltonian can then be found in a perturbative expan¬ 
sion iJeff = -I- iLgg -I- ... in powers of using 

different methods [31-34]. The lowest-order term 
corresponds to the average or static Fourier component 
of the Hamiltonian. Higher-order terms, on the other 
hand, depend on the particular choice of gauge Uf{0)- 
For convenience, we choose the gauge [35] that leads to 
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with the Fourier components ^ 

We consider spinless, non-interacting particles on a 
shaken hexagonal lattice described by the Hamiltonian 
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= E Cr. + Cr. + H.c.)), 
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with the vector creation and annihilation operators c\. = 
and Cr- = {cA,i,CB,i)"’" satisfying the usual 
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(anti)commutator relations, and the time-dependent ma¬ 
trices 
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with energy offset A. The summation in Eq. (3) is per¬ 
formed over the positions of all unit cells of the hexag¬ 
onal lattice, which we consider to be infinite or with pe¬ 
riodic boundary conditions. The vectors bi = a(-\/3, 0) 
and b 2 = |(—3) correspond to two primitive vec¬ 
tors that span the underlying triangular Bravais lat¬ 
tice, with the distance a between nearest-neighbor (NN) 
sites. The time-dependent rates that characterize the 
NN tunneling are given by 9ak{t) = jo with 

Xk{t) = fgdr F(t) ■ afe - X fj'f* dr F(r) • in 
terms of the driving force F(t) and the real tunneling 
amplitudes jo of the undriven system. In general, the 
time-dependent tunneling rates depend on the direction 
of tunneling defined through the vectors that connect 
neighboring sites ai = K-s/S, 1 ), a 2 = |(—-x/S, 1 ) and 
as = —ai — a2. 

As theoretically expected and experimentally con¬ 
firmed [ 8 ], the dynamics of the system in a fast driving 
regime can be captured very well by the truncated effec¬ 
tive Hamiltonian 

i?dh = i?eff^+i?eff- (6) 


Since the leading-order effective Hamiltonian is given by 
the average of ff(t) in Eq. (3), contains the same 
tunneling processes as the undriven system: the on-site 
energies ±A remain invariant and the effective tunnel¬ 
ing rates become the directionality-dependent quantities 

dL = tIo 

The first-order term given by Eq. (2) reads 

Cr.+H.C., (7) 

i k—0 

where bo = 0, bs = —bi — b 2 . The effective matrices 
= diag(T/c,-Tfc) with tq = I]i=i '^i = 

?ii(a 2 ,-a 3 ) T 2 = w(a 3 ,-ai) and T 3 = ?ii(ai,-a 2 ) are 
defined in terms of 


00 . 
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FIG. 1. Phase diagram of the isotropic effective Hamilto¬ 
nian in Eq. (9) (black lines and dark colors) overlapped with 
the phase diagram of the Haldane model [36] (gray lines and 
lighter colors), giving the Chern number C of the lowest ern- 
ergy band as a function of the phase (f) and ratio A/j 2 . Or¬ 
ange represents C = — 1, blue C = 1 and white C = Q. For 
(j> = ±7r/2 the Chern number of the two Hamiltonians coincide 
independently of A/jh. 


those of the Haldane model [7], where the two tunneling 
rates are complex conjugated with respect to each other. 
Only for purely imaginary rates = — r^, k = 1,2,3, 
do the NNN rates of the two models coincide. Conse¬ 
quently, it is fundamentally impossible to implement the 
full topological diagram of the Haldane model via lattice 
shaking without non-vanishing real NNN tunneling rates 
in the undriven system. 

Despite the differences between the Haldane model 
Hamiltonian and the effective Hamiltonian in Eq. ( 6 ), 
the two models share similar topological properties. As 
we shall later demonstrate, it is possible to find a driving 
force yielding isotropic tunneling rates, namely 5 °^ = ji 
and Tfc = for all directions k = 1 , 2 ,3, where ji and 
j 2 are positive real numbers and (j is defined in the inter¬ 
val (—TT, tt]. The effective Hamiltonian can then be writ¬ 
ten in quasimomentum space as H^h = 
where are the vector momentum creation and anni¬ 
hilation operators and 


3 

i7(k) = ^h,(k)(T, (9) 

is defined in terms of the Pauli matrices Ui and 


with the Fourier components ? /o^ 

The rates Tk, k = 1,2,3, describe effective next-nearest- 
neighbor (NNN) tunneling processes that result from a 
virtual tunneling process over a neighboring site. The 
relative sign in between the different rates and 
—Tk is a fundamental symmetry that is independent of 
the specific driving force F{t). 

Due to this symmetry, the emergent NNN tunneling 
rates discussed above are, in general, not equivalent to 


/ii(k) = ji (1 -I- cos(k • bi) + cos(k • b 2 )), (10) 

/i 2 (k) = ji (sm(k • bi) - sin(k • b 2 )), ( 11 ) 

3 

/i3(k) = A-k 2j2 E 

i=l 

The topological diagram of this model, displaying the 
Chern number [5] as a functions of the Hamiltonian pa¬ 
rameters A/j 2 and (j, can be readily calculated [37] and 
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it is shown in Fig. 1. The transition between differ¬ 
ent topological phases - indicated with a solid black 
line - corresponds to parameters of the Hamiltonian for 
which the gap between the two energy bands e±(k) = 
^^h‘l + h 2 + hi closes. For comparison, we also display 
the analogous topological diagram of the isotropic Hal¬ 
dane model Hamiltonian [36]. Only for (p = ± 7 r /2 do the 
two Hamiltonians coincide, consistently with the diagram 
in Fig. 1. 

In order to assess to what extend the entire parameter 
regime of the topological diagram can be realistically ex¬ 
plored, it is necessary to correctly identify driving forces 
F(t) that lead to isotropic effective tunneling rates with 
controllable amplitudes and phase p. Since the topolog¬ 
ical energy bands emerge as a consequence of the inter¬ 
play between the NN and NNN tunneling processes, it is 
important that the relative effect of the NNN tunneling 
with respect to NN tunneling, given by the ratio jij is 
sufficiently large. Nevertheless, these two tunneling pro¬ 
cesses are of different orders of magnitude, since ji ~ jo 
andj 2 ^ As the ratio jd/ji is proportional to jo/w, 
it could easily be increased through a decrease of the driv¬ 
ing frequency w. This, however, could compromise the 
validity of the high frequency expansion of the effective 
Hamiltonian. For this reason, we consider a small fixed 
ratio jo/w to be determined according to the experimen¬ 
tal setup and aim at finding a driving force with a set of 
parameters p that maximize the proportionality factor 
between j 2 /ji and jo/w. Since the amplitude ji is 
directly related to 0 , we introduce p = (jtg as a constraint 
for the maximization, where tjtg is the desired phase that 
we target. Additionally, ji should be sufhciently large 
with respect to the bare tunneling rate jo in order to 
avoid that the effective tunneling processes appear at the 
expense of slowing down the dynamics as compared to the 
undriven system. We therefore introduce the additional 
constraint ji/jo > rth, where the threshold value rth can 
be chosen from the interval 0 < rth < 1 - 

We thus aim at finding a driving force targeting: 

(i) Isotropy 5 °^ = ji and = jse*'^, fc = 1, 2, 3. 

(ii) Controlability of the phase p. 

(hi) Enhancement of the NNN tunneling rates through 
the maximization 


R{4>tg,rth) = jmax^— ^ > rth ; = <jtg| (13) 

I P JiJo Jo J 

performed over a set of free driving parameters p. 

For a monochromatic driving force [ 8 , 38-40], the effec¬ 
tive NNN tunneling rates become purely imaginary and, 
thus, only the points 4> = ± 7 r /2 in Fig. 1 can be accessed. 
This strong limitation can be overcome by specifically 
designing the driving pulse to satisfy our requirements. 
Despite the highly non-linear dependence that the tun¬ 
neling rates have on the driving force, we find analytic 
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FIG. 2. Phase <f) of the complex effective next-nearest- 
neighbor tunneling rates as a function of AxjLo and Aiju) for 
a force F+(f) with N = 2 and 82 = 7 r/ 2 . All phases can be 
explored with a suitable choice of Ai and A 2 . The dashed and 
solid white lines indicate the contour lines for jr/jo = 0.5 and 
ji/jo = 0.25 and limit the region of accessible phases given a 
constraint with rth = 0.5 or rth = 0.25, respectively. 



expressions for and rj, in terms of driving parameters 
by using Wdimensional Bessel functions [41]. 

This allows us to identify the structure of two driving 
forces F 4 .(t) and F_(t) that lead to isotropic tunneling 
rates independently of their free parameters [41] . Conse¬ 
quently, the requirement (i) above is automatically satis¬ 
fied [42] , which allows the examination of the subsequent 
target properties. The general form of the driving forces 
F±(t) containing N different frequency harmonics reads 

N 

F±(t) = X! (cos(w„t - (5„)ei -h cos{ujnt - 6 ^) 62 ) 

n—1 

(14) 

with the two perpendicular vectors ei = (ai — a 2)/-\/3 
and 02 = —aa, the relative phases 5 ^ = ± (—l)"7r/2 

and the frequencies ujm = \ (6m — (—1)"* — 3)cj, which 
parametrize all positive integer multiples of w except 
those that are multiples of 3cij. Because the overall phase 
of the driving force is irrelevant in the fast-driving regime, 
we choose = 0 in the following. The remaining 2N — 1 
driving parameters comprise the set p and need to be 
chosen so that the requirements (ii) and (iii) are satis¬ 
fied. 

In order to ease an experimental implementation, we 
consider the simplest polychromatic force with N = 2, 
which contains three driving parameters: Ai, Ai and Si. 
We analytically find that if <52 = 0, ±7r, the real part of 
Tfe again vanishes, similarly to the monochromatic case. 
However, for 5i ^ 0, ± 7 r and appropriate choice of driv¬ 
ing amplitudes the entire range of phases 4> G (—tt, tt] 
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FIG. 3. Plot of the real and imaginary part of R{4>tg, 
as a function of a discrete set of target phases 0tg and for two 
different values of rth- In polar coordinates, the radius and 
argument of each data point coincide with R{(j>tg,rth) and 
4>tg, respectively. Dots correspond to rth = 0.25, and crosses 
to rth = 0.5. The results in blue have been obtained with a 
driving force F+(t) and the results in orange with F_(t). 


can be realized, satisfying thus the requirement (ii), as 
illustrated in Fig. 2. 

The maximization in the point (iii) leads to a signifi¬ 
cant enhancement of the effective NNN tunneling for any 
desired phase, but considerably different results are ob¬ 
tained depending on the phase (ptg and threshold rate rth 
that we target. In order to discuss this behavior, we plot 
in Fig. 3 the real and imaginary part of 
for two different values of rth and for a discrete set of 
angles 0tg € (~'^i7r]. Each data point is obtained with 
a numerical constrained optimization over the set of free 
parameters p = {Ai/uj,A2/uj,S2)- Overall, we observe 
two main features. 

First, for a given threshold value rth, the largest val¬ 
ues of i?(0tg,Dh) are obtained for a range of phases 
close to ±7r/2. The lowest values correspond to phases 
0 and ±7r, for which effective Hamiltonian Ff^h is time- 
reversal invariant. This indicates that experimentally it 
is easier to access the areas of the topological diagram 
in Eq. (1) that are near (p = F-k Noteworthy, we 
find that for c^tg = ±7r/2 the optimal solution of the two- 
frequency pulse reduces to a monochromatic force, i.e. 
A 2 = 0, independently of rth. Nonetheless, the maximum 
of i?((?!)tg, rth) does not always correspond to ptg = ±7r/2 


for a fixed rth, as can be seen in the results for rth = 0.25 
in Fig. 3. 

Second, the lower the threshold value rth is for a fixed 
target phase ptg, the larger R{4>tg, Dh) can be, as a larger 
region in the parameter space, given by the driving pa¬ 
rameters, can be accessed (see contour lines in Fig. 2). 
Thus, there is a trade-off between lower values of the 
threshold rth for the ratio ji/jo and higher relative en¬ 
hancement 4^ of NNN tunneling. It is thus advisable to 

Ji Jo 

choose a small value of rth, provided that it is sufficiently 
large so that tunneling is dominant in the dynamics of 
the system and processes like interaction or heating can 
be neglected on the time-scale on which tunneling occurs. 

As a result of the maximization in Eq. (13), the set 
of optimal driving parameters that yield the maximum 
enhancement of the NNN tunneling are also obtained. 
We find that for all the results shown in Eig. 3, the 
corresponding driving amplitudes are of a similar order of 
magnitude as the driving frequency, specifically |Ai/a;| < 
3.5, * = 1,2. Additionally, we observe a discontinuous 
behavior of the optimal driving parameters as a function 
of 4 >tg, which leads to a discontinuity in R{(j)tg,rth), see 
rth = 0.25 results in Fig. 3. This can be understood 
in terms of the constraint jxjjo > Dh, which restricts 
the parameter space to disjoint regions, as can be seen 
in the contour line ji/jo = 0.25 in Fig. 2. Depending on 
the targeted phase, the driving parameters might change 
from one region to another, yielding a discontinuity in 
the driving parameters and in the corresponding value of 
R{(ptg, Dh)- 

Increasing the number of frequency harmonics of the 
force F±(t), more driving parameters become accessible. 
These additional degrees of freedom can be exploited e.g. 
to further increase the maximum value of i?(0tg, "Cth) tar¬ 
geting specific phases, as we demonstrate for a driving 
force with N = 3 [41]. 

Our results show that even a very low number of fre¬ 
quency components is sufficient to completely outperform 
the monochromatic driving and yield significant NNN 
tunneling rates for any phase (p. This exemplifies the 
fact that the usually-considered monochromatic driving 
can strongly limit the accessible effective dynamics and 
that only suitably chosen driving protocols can enable 
the exploration of the entire accessible dynamics. Fi¬ 
nally, we believe that the possibility to experimentally 
implement the model in Eq. (9) in deep optical lattices 
naturally warrants further theoretical investigations, e.g. 
on the emergence of fractional quantum Hall effect [43] 
when interactions are taken into account. 

We are indebted to Julian Struck, Juliette Si- 
monet, Gregor Jotzu, Andreas Mielke and Elizabeth 
von Hauff for stimulating discussions or feedback on 
the manuscript. Financial support by the European 
Research Council within the project ODYCQUENT is 
gratefully acknowledged. 
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This Supplemental Material is divided into three sec¬ 
tions. In Sec. I, we introduce the mathematical tools 
necessary to derive analytical expressions for the tunnel¬ 
ing rates of the effective Hamiltonian of shaken lattices. 
In Sec. II, we demonstrate that the driving forces in¬ 
troduced in Eq. (14) lead to isotropic tunneling rates. 
Finally, in Sec. Ill we provide the optimal results corre¬ 
sponding to a force with three Fourier components. 

I. MULTIDIMENSIONAL BESSEL FUNCTIONS 

The effective tunneling rates of the Hamiltonian 
in Eq. (6) are given in terms of the Fourier components 
g2 ., which are defined through 

OO 

n— — oo 

with a T-periodic function Xkit)y as described in the main 
text. Motivated by the need of analytical expressions to 
target the properties described in the points (i)-(iii), we 
introduce multidimensional Bessel functions as a general¬ 
ization of ordinary and two-dimensional Bessel functions 
[ 11 - 

Ordinary n-th order Bessel functions can be de¬ 

fined as Fourier components of the generating function 

OO 

g*asin(h) ^ (S2) 

n——oo 

Here, we generalize the Bessel functions introduced above 
and consider A^-dimensional Bessel functions Jnic, d) de¬ 
fined by the generating function 

OO 

= Y J'„^(c;d)e“^, (S3) 

n——oo 

with the real vectors c = (ci,...,CAr) and d = 

(di,..., dAf), and the integer vector z = { zi ,..., zn )- 
The Fourier transform of Eq. (S3) leads to the integral 


representation 

1 

J^{c;d) = — / (S4) 

27r Jo 

Inserting Eq. (S2) into Eq. (S4) permits to expand N- 
dimensional Bessel functions in terms of ordinary one¬ 
dimensional Bessel functions 

J„^(c;d) .Jx„(c^)e-*‘'" 

X 

— / (S5) 

2'7r Jq 

= ^^xi(ci).J.^(cjv)e-''^"dz.x,„(S6) 

X 

= Y.JsAci) .J.„(cw)e-*'^•^ (S7) 

S 

where x = {xi ,..., xn) denotes all possible integer vec¬ 
tors and s = (si,..., sn) all integer solutions of the Dio- 
phantine (i.e. integer) equation 

z • s = ziSi -I— • -I- znSn = n. (S8) 

A linear Diophantine equation of the form in Eq. (S8) 
can be solved if and only if the greatest common divi¬ 
sor gcd( 2 ;i ,... ,zn) divides n [2]. Since the driving force 
will always contain the fundamental driving frequency, 
at least one of the elements of z equals 1, which implies 
that Eq. (S8) can be trivially solved for all n. Without 
loss of generality we order the elements such that Zi = 1. 
The solution of Eq. (S8) then reads 

N 

(S9) 

i^2 

where S 2 ,..., sat are N — 1 integer free parameters that 
characterize the general solution. Consequently, since 
Zi = 1, A-dimensional Bessel functions can be expressed 
as 


CJO 

J:(c;d)= ^ .(SIO) 


A numerical evaluation of N-dimensional Bessel func- tions can be implemented by truncating the infinite sums 
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in Eq. (SIO). Specifically, in the numeral implementa¬ 
tion of the maximization described in the main article we 
examine amplitudes with magnitudes \ci\ < 5, for which 
a truncation of |si|, z = 1, • • • , up to ~ 7 is sufficient. 
We find that, over a wide range of parameters, this im¬ 
plementation is computationally more efhcient than other 
numerical approaches such as a Taylor series truncation 
of the generating function in Eq. (S3). Nonetheless, the 
dimension of the Bessel fuirction rapidly becomes a lim¬ 
iting factor since the number of terms to evaluate grows 
with N. For this reason, it is convenient to work with the 
lowest dimension possible, which is given by the number 
of elements in z that are different. If a subset z' of z 
contains the same elements, i.e. z' = (z,..., z), then the 
polar parametrization 

re—= ^c'e-*< (Sll) 

i 

allows to express 

^ c'sin(zr — d') = r sin(zT — a), (S12) 

i 

where r > 0 and c', d' are the elements of c and d as¬ 
sociated to z'. In this manner, the dimension of z can 
be reduced by dim(z') — I. This process can be repeated 
until all the elements of z are different. 


II. ISOTROPIC TUNNELING RATES 

In this section we will demonstrate how the driving 
forces F±(t) introduced in Eq. (14) lead to isotropic 
tunneling rates and defined after Eq. ( 6 ) and (7), 
respectively. For this purpose, we shall first find ana¬ 
lytical expressions in terms of multidimensional Bessel 
functions for the Fourier components of the time- 
dependent tunneling rates in Eq. ( 8 ). 

We start by considering the force F+(t). The quantity 
Xk{t) introduced after Eq. (5) then reads 


Moreover, from Eq. (S15) we find that the directionality- 
dependent phases can be expressed as 

dt = S„, + {-ir^k, (S16) 

where (j)k denotes the angle of the vector with respect 
to ei and they read 4>i = 7 r/ 6 , 4>2 = Stt/G and ^3 = 37 r /2 
with the convention introduced in the main text after Eq. 

( 5 ). ... 

Using this representation, the number of terms in Eq. 
(S13) is reduced from 2N to N such that 

N 

Xk{t) = X! sin(w„t - d^). (S17) 

m—1 

With isotropy in the undriven system, i.e. = jg, 

fc = 1, 2,3, the Fourier components can be expressed 
in terms of A-dimensional Bessel functions as 


g:,=joj:ic-,d^), (SI 8 ) 

where c = (ci • • • ,CAr), = (d^ • • • ,d%) and z = 
(zi, • • • , zjv) of the Bessel functions, with z^ = 

The directionality dependence of the Fourier compo- 

- j fc 

nents appears only in the exponential * of the 

Bessel function, see Eq. (SIO). In fact, due to Eq. 
(S16), the only possible directionality dependence origi¬ 
nates from = {—4>k, ■ • ■ , (—l)'^</>fc) through 


N 


N 


e *P ■"* = exp z n - ^ ZmSm (t>k - z ^(-l)™^fcS,i 


N 


= exp z 0 fc n - ^ (z„ -b (-l)'")s 


mil? 


(S19) 


m=2 


with the integer vector 


N 


— In ^ ‘ ‘ * 5 

\ m=2 J 


(S20) 


From Eq. (S19), the isotropy of the nearest-neighbor 
tunneling rates directly follows. For rz = 0, Eq. (S19) 
reduces to 


N 

Xkit) = ^ - Sm)a-k ■ ei 

m=l 

+ sm{uJmt - 5^)ak ■ 62 ). (S13) 

Following the same arguments as in the previous section, 
we introduce the polar representation 

= — (e“*'^"*afc • ei -b e“*^-afe • 62 ) (S14) 

= f^e-*^™(a,. ■ ei - z(-l)"'afc • e2).(S15) 

Importantly, since the vectors ai, a 2 and as have the 
same amplitude and because the chosen relative phases 
Sm'^ satisfy — dm = (—l)™7r/2, the quantities Cm do 
not depend on the direction specified by the index k. 


-ip^-s ( -OJ. (^m + (-l)™) ^ ZC01^ 

e P = exp -z3d)fe ^ ---Sm (821) 

\ m=2 / 

Because of our choice of ujm in Eq. (14), the quantity 
(Zm + (—I)"*) appearing in Eq. (S19) is always a multi¬ 
ple of 3 and thus (zm + (—1)™)/3 is an integer number. 
Consequently, since the phases (j)k satisfy 

^-i34>iq _ g-i3b>29 _ g-j3<b39 (S22) 

for arbitrary integers q, we obtain that Jo(c;d^) and 
hence are independent of the direction k. 

Finally, we will demonstrate the isotropy of the next- 
nearest-neighbor tunneling rates Tfc by demonstrating 
that 

_ g-ip^ Sg-ip^-s' 

• 3 - 1 / 


(S23) 
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with the integer vectors s in Eq. (S20) and 


which implies 


N 


s' = -n - ^ s'r 


■’2) ) I ) 


(S24) 


J:(c; di) d^) = J^{c; d^) JZJc; d^) = J'„(c; d^) JL„(c; d^), 


(S25) 


and hence ti = T 2 = T 3 . 

In order to prove Eq. (S23), expressions for e“*P * are necessary. They can be obtained with Eq. (S19) and 

read 


N 


o-*P Sp-ip^S 


= exp -i3(^i ^ 


{z^ + i-ir) 


N 


exp ^ 


(Zm + { !)'") „/ j 


(S26) 


m—2 


m—2 


Eq. (S23) is then derived from Eq. (S26) by using 
Eq. (S22) and the relation = 

gin( 03 - 0 i)^ 

The same result can by readily derived for F_(t) by 
noting that F_(t) is obtained from F_(.(t) with the sub¬ 
stitution ui —>■ —u! and Sn —?> —Sn- 


III. OPTIMAL RESULTS WITH THREE 
FOURIER COMPONENTS 


In the main Letter, we have shown how a suitably 
chosen driving force F ± (t) with two Fourier components 
leads to an optimized finite ratio for any desired 

phase 4>. Here, we demonstrate an enhancement of such 
ratio by increasing the number of frequencies to A = 3 
and properly choosing the extra driving parameters. 

A three-frequency driving force contains a set of driv¬ 
ing parameters p = (Ai /w, A 2 /w, A 3 /w, ^2 , 1^3 ), where A 3 
and ^3 correspond to the amplitude and phase associated 
to the frequency component = Aio. Analogously to the 
two-frequency results, we obtain the optimal choice of p 
with the constrained maximization in Eq. (13) [3], with 
the target functional and constraints expressed analyti¬ 
cally in terms of 3-dimensional Bessel functions. In Fig. 
1 , we display optimal results obtained and compare them 
with the N = 2 results. We observe that, for the entire 
range of target phases (ptg except for (j)tg = ±7r/2, op¬ 
timally chosen forces with N = 3 lead to larger values 
of i?(0tgUth) than for bichromatic driving. Moreover, 
similarly to the N = 2 case, we find that the optimal 
driving force targeting (j)tg = ±7r/2 corresponds to the 
monochromatic force, i.e. A 2 = A 3 = 0. 

Remarkably, the enhancement of the effective next- 
nearest-neighbor tunneling as compared to the nearest- 


neighbor tunneling, that is, the enhancement of effects 
of order ~ as compared to effects of order ~ 1, re¬ 
lies on the use of a higher frequency component. This 
strongly suggests that the use of driving forces F± (t) with 
even higher frequency harmonics, that is with N > 3, 
and optimally chosen driving parameters could still fur¬ 
ther improve on the value of R{(j)tg,rth)- Nonetheless, 
our results show that already a very low number of fre¬ 
quency components is sufficient to significantly outper¬ 
form the monochromatic driving and yield significant 
next-nearest-neighbor tunneling rates for any phase <f>. 
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FIG. 1. Plot of R(0tg,rth = 0.5) calculated numerically as a 
function of a discrete set of target phases <f)tg = <j)n = n-K j40 
indexed by n = 1, • ■ • , 20. Thus, the label n — 0 and n = 20 
correspond to 0tg = 0 and </>tg = '7r/2, respectively. In green, 
results corresponding to a driving force F+(t) with N — 3. 
For comparison, with blue crosses we display again the results 
for N = 2. 
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demic Press, 1969). 

[3] In practice, we do not exactly require constraint (j) = (fitg 


but rather \(f> — (f>tg\ < 5(f), where 5(fi = 0.01 is the error 
that we admit. 



